arXiv:1507.00596v2 [math.NA] 9 Dec 2016 


A fast and well-conditioned spectral method for singular integral equations 


Richard Mikael Slevinsky^ *, Sheehan Olver'’ 

^Department of Mathematics, University of Manitoba, Winnipeg, Canada 
^School of Mathematics and Statistics, The University of Sydney, Sydney, Australia. 


Abstract 

We develop a spectral method for solving univariate singular integral equations over unions of intervals by utilizing 
Chebyshev and ultraspherical polynomials to reformulate the equations as almost-banded infinite-dimensional sys¬ 
tems. This is accomplished by utilizing low rank approximations for sparse representations of the bivariate kernels. 
The resulting system can be solved in OirrP'n) operations using an adaptive QR factorization, where m is the bandwidth 
and n is the optimal number of unknowns needed to resolve the true solution. The complexity is reduced to 0{mn) 
operations by pre-caching the QR factorization when the same operator is used for multiple right-hand sides. Stability 
is proved by showing that the resulting linear operator can be diagonally preconditioned to be a compact perturbation 
of the identity. Applications considered include the Faraday cage, and acoustic scattering for the Helmholtz and grav¬ 
ity Helmholtz equations, including spectrally accurate numerical evaluation of the far- and near-held solution. The 
Julia software package SingularIntegralEquations. jl implements our method with a convenient, user-friendly 
interface. 

Keywords: Spectral method, ultraspherical polynomials, singular integral equations 
2010 MSC: 65N35, 65R20, 33C45, 31A10. 


1. Introduction 

Singular integral equations are prevalent in the study of fracture mechanics fTl, acoustic scattering problems ^ 
EElElIll, Stokes How Q, Riemann-Hilbert problems il, and beam physics HSlIIl- We develop a fast and stable 
algorithm for the solution of univariate singular integral equations of general form IfTTII 

K{x,y)u(y)dy = fix), for .r e T, Su - c, (1) 

where K{x, y) is singular along the line y - x, the x in the integral sign denotes either the Cauchy principal value or 
the Hadamard Hnite-part, F is a union of bounded smooth open arcs in K^, and S is a list of functionals. To be precise, 
we consider the prototypical singular integral equations on [- 1 , 1 ] given by: 

“"T + = /(v), for xe[-l,l], 

y-x j 

where Ai,..., A 4 are continuous bivariate kernels. 

In this work, we use several remarkable properties of Chebyshev polynomials including their spectral convergence, 
explicit formulae for their Hilbert and Cauchy transforms, and low rank bivariate approximations to construct a fast 
and well-conditioned spectral method for solving univariate singular integral equations. Chebyshev and ultraspherical 
polynomials are utilized to convert singular integral operators into numerically banded inhnite-dimensional operators. 
To represent bivariate kernels, we use the low rank approximations of IfT^ . where expansions in Chebyshev polyno¬ 
mials are constructed via sums of outer products of univariate Chebyshev expansions. The minimal solution to the 
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recurrence relation is automatically revealed by the adaptive QR factorization of ifTJII . Diagonal right preconditioners 
are derived for integral equations encoding Dirichlet and Neumann boundary conditions such that the preconditioned 
operators are compact perturbations of the identity. 

The inspiration behind the proposed numerical method is the ultraspherical spectral method for solving ordi¬ 
nary differential equations E], where ordinary differential equations are converted to infinite-dimensional almost 
banded linear systems (an almost banded operator is a banded operator apart from a finite number of dense rows). 

These systems can be solved in infinite-dimensions, i.e., without truncating the operators m, as implemented in 
ApproxFun.jl ca in the Julia programming language 016111711 . The Julia software package SingularIntegralEquations . j 1 E 
implements our method with a convenient, user-friendly interface. As an extension of this framework for infinite¬ 
dimensional linear algebra, mixed equations involving derivatives and singular integral operators can be solved in a 
unified way. 

Several classical numerical methods exist for singular Fredholm integral equations of the first kind. These include: 
the Nystrom method imiaoiEii, whereby integral operators are approximated by quadrature rules; the collocation 
method EilEa, where approximate solutions in a hnite-dimensional subspace are required to satisfy the integral 
equation at a hnite number of collocation points; and the Galerkin method ll24ll^ . where the approximate solution is 
sought from an orthogonal subspace and is minimal in the energy norm. The use of hybrid Gauss-trapezoidal quadra¬ 
ture rules ll26ll27ll28ll29l can signihcantly increase the convergence rates when treating weakly singular kernels. 

Numerous methods have exploited the underlying structure of the linear systems arising from discretizing inte¬ 
gral equations. The most celebrated of these is the Fast Multipole Method of Greengard and Rokhlin ll^ . Other 
characterizations in terms of semi-separability or other hierarchies have also gained prominence ||3T] |32] |33l . Ex¬ 
ploiting the matrix structure allows for fast matrix-vector products, which then allows for Krylov subspace meth¬ 
ods ll^ to be extremely competitive. For scattering of the Helmholtz equation in very special geometries, hybrid 
numerical-asymptotic methods have been derived for frequency-independent solutions to the Dirichlet and Neumann 
problems ll4l [35l 1^ 1^. 

Previous works on Chebyshev-based methods for singular integral equations include Frenkel llTTll . which derives 
recurrence relations for the Chebyshev expansion of a singular integral equation after expanding the bivariate kernel 
in a basis of Chebyshev polynomials of the hrst kind in both variables, and Chan et al. Il^[^ in fracture mechanics, 
among others. A similar analysis in Il40ll is used for hypersingular integrodifferential equations by expanding the 
bivariate kernel in a basis of Chebyshev polynomials of the second kind. This paper is an extension of these ideas 
with essential practical numerical considerations. 

Remarks. 

1. Combined with fast multiplication of Chebyshev series, our method is suitable for use in iterative Krylov sub¬ 
space methods. 

2. There is a great diversity of integral equation formulations. The choice of formulation depends on many prop¬ 
erties, including for example, whether the boundary is open or closed and whether there are resonances. Most 
equations involve operators that contain manipulations of the fundamental solution, which would still satisfy the 
requirements of our method. However, we focus on the direct integral equations to retain a simple exposition. 


2. Boundary integral equations in two dimensions 


In two dimensions, let x = (xi,X 2 ) and y = (ji,y 2 )- Positive definite second-order linear elliptic partial differential 
operators (PDOs) with variable coefficients are always reducible to the following canonical form IHTI : 


A 

L|m) — Au + a- -h b- -h cu. 

dx\ dx 2 


( 2 ) 


Let 0(x, y) denote the positive definiXe fundamental solution of (|^ satisfying the formal partial differential equation 
(PDE) 


LaI®) = -b(x - y). 


( 3 ) 


where 6 is the two-dimensional Dirac delta distribution and the subscript indicates that L is acting in the x variable. 
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2.1. Exterior scattering problems 

Let r be a union of disjoint bounded smooth open arcs in and let Q := \ L. 

Let m(x) dy be the standard Fourier transform in Then for i e M, defines the Bessel 

potential space as the set of notmed tempered distributions equipped with: 

II«IIto= f (l + |x|"yi«(x)pdx. (4) 

JR2 

For bounded Q with non-empty interior, the space FF'(Q) := jwln : u e Furthermore, let denote the 

spaced of locally integrable functions in FF'(Q) and lastly, let H^{Q.) \u e : suppM c flj. 

Definition 1. For x e Q, let .Sr and Dr define the single- and double-layer potentials: 

.SrM(x) = 0(x, y)M(y) dF(y) : H ^ (F) ^ hIJQ), (5) 

^)rM(x) = £ ^ 

Definition 2 (Radiation condition at infinity ||42|| ). We say that u e satisfies the radiation condition at in¬ 

finity if: 

lim (.S|y|=p [duldn] (x) - D\y\=p [u] (x)| = 0, for x 6 Q. (7) 

For solutions of the homogeneous equation L{m) = 0 satisfying the radiation condition at infinity. Green’s repre¬ 
sentation theorem allows for the determination of the exterior solutions given data on the boundary F: 

m(x) = -.Sr [dujdn] (x) + Dr [u] (x), for x 6 Q. (8) 


Here, [m] denotes the jump in u along F and [dujdn] the jump in its normal derivative. These are formally defined 
by the Dirichlet trace and conormal derivative fTS), or in the case of the Laplace equation, simply as the difference 
between the limiting values on F as we approach from the left and the right. This identity can be interpreted as 
representing u in terms of the potential of a distribution of poles on F through the single-layer and normal dipoles on 
F through the double layer. With either Dirichlet or Neumann boundary conditions, we restrict (|^ to the boundary 
and solve for the unknown boundary value. Once both quantities on the boundary are determined, the solution to the 
exterior problem is readily available in integral form. 

Dirichlet Problem Given an incident wave m'(x) e H^i^) satisfying L{m') = 0, find m^(x) e satisfying 

L{m') = 0, the radiation condition at infinity, and 

m'(x) h- u \ x ) = 0, for X e F. (9) 


Neumann Problem Given 


du'ix) 


dn(x) 

radiation condition at infinity, and 


e H 2 (F) satisfying L{m') = 0, find m‘'(x) e satisfying L{m^) = 0, the 


d 

dn(x) 


(m'(x) 


-H M^(X)) = 0, 


for X e F. 


( 10 ) 


For the case of the Laplace and Helmholtz equations, the Dirichlet problem is originally formulated in ll43l Eqs. 
(1.1) & (1.2), (1.6) & (1.7)]. Similarly, the Neumann problem is originally formulated in P4] Eqs. (1.1) & (1.2)]. 

Dirichlet Solution The Dirichlet problem is solved by (|^ where [m^] = 0, and the scattered solution is represented 
everywhere by the single-layer potential. The density [duydn] 6 in ^ satisfies: 


X 


‘l>(x,y) 


du” 

dn 


dr(y) = u‘(x), for x e F. 


( 11 ) 
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Neumann Solution The Neumann problem is solved by ([^ where {du^ ldn\ - 0, and the scattered solution is 
represented everywhere by the double-layer potential. The density [m^] e (T) in ([^ satisfies: 


d r d(I>(x, y) 
dn(x) Jr dn(y) 


K]dr(y) = - 


du‘(x) 
dn(x) ’ 


for X e r. 


( 12 ) 


For the case of the Laplace and Helmholtz equations, the Dirichlet solution is originally proved in Il43l Theorems 
1.4 & 1.7]. Similarly, the Neumann solution is originally proved in Il44l Theorem 1.3]. Furthermore, by appealing 
to the theory of Mellin transforms, inverse square root singular behaviour is derived for the open ends of F in the 
Dirichlet problem ll43l Theorem 2.3], and square root singular behaviour is derived for the open ends of F in the 
Neumann problem ll44l Theorem 1.8]. 

For elliptic PDOs with variable coefficients, we use the solutions provided by the singular integral equations ( [TT] i 
and ([T2]i, and while the scope of this paper is numerical, we conjecture that they hold more generally. 


2.2. Riemann functions 

In addition to the PDO in (|^, consider its adjoint: 


T*r 1 A , 

L 1 y) = Av-^--h cv. 


dxy dX2 


With the change to complex characteristic variables: 

z — xi + ix 2 , ^ - xi - ix 2 . 


zo = yi + m, ^0 = yi - m, 


(13) 


(14) 


L and L* take the form: 


where: 


L{t/] = 
L*{y] = 



d^U 

dU 

dU 


A— +B — +CU. 

dzd^ 

dz 

d( 

d^V 

diAV) 

, cv 

dzd^ 

dz 


/z + ^ 


+ 

+ 

1 

\ 2 ’ 

2i / 

V 2 2i 

(z + ^ 


1 

+ 

1 

\ 2 ’ 

2i / 

V 2 2i 

z + ^ z 



2 ’ 

2i /■ 



) 

) 


(15) 

(16) 

(17) 

(18) 
(19) 


Theorem 3 (Vekua and Garabedian ||4T]|45l)- For analytic functions ([ni'-'lIEI. there exist analytic functions 9\{z, f, Zo, ^o) 
and go(z, Zo, ^o) such that: 


4>(z, ^,zo,fo) 


1 

An 


(H(z, f, zo, ^o) log[(z - zo)(f - ^o)] + go(z, f, zo, ^o). 


( 20 ) 


where L{<1)] = 0 in (z, f) and L*{<E)] = 0 in (zo, ^o) so long as z ^ Zo and ^ + fo- In 93 is the Riemann function of 

the operator L satisfying: 


L*{93] = 0, 
93(zo,^,zo,Co) = exp 

93(z,^o,2o,^o) = exp 



A(zo, t) drj 

B(f,^o)df|. 


and 


( 21 ) 

( 22 ) 

(23) 
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Remarks. It is straightforward to reformulate (|2T]i-(|27| to the following integral equation: 


f 




ff 

■Jzo J(o 


- f A(z,T)D^(z,T,zo,^o)dT 
J(o 

C{t,T)Dl(t,T,zo,^o)dTdt = 1. 


( 24 ) 


Returning to the original coordinates x and y, fundamental solutions for elliptic PDOs with analytic coefficients 
can be written as: 

(l)(x, y) = A(x, y) log |x - y| + B(x, y), (25) 

where A and B are both analytic functions of x and y and where A(x,y) = -^Dl(z,^,zo,^o) implying A(x,x) = 
-(27r)“'. If, furthermore, the PDO is formally self-adjoint, then A and B are also symmetric functions of x and y. 


3. Practical approximation theory 

Chebyshev approximation theory is a very rich subject that has seen numerous exceptional contributions: see ll46l 
@7]51] and the references therein. In this section, we describe some approximation spaces for one-dimensional inter¬ 
vals and two-dimensional squares. For every approximation space, one may consider the interpolants, which are equal 
to the function at a set of interpolation points, and the projections, which are truncations of the function’s expansion. 
Unless an extraordinary amount of analytic information is known about a function, interpolants are generally easier 
to construct. 

We consider an approximation space practical if there is a fast way to transform the interpolation condition into 
approximate projections. While a few methods exist to create fast transforms, all the practical approximation spaces 
we consider resort to some variation of the fast Fourier transform (FFT) ll49l l50l to reduce OirP') complexity to 
0{n\ogn). Other properties which make an approximation space practical are: 0{n) evaluation; a low Lebesgue 
constant; absolute, uniform, and geometric convergence with analyticity; and, easy manipulation for the development 
of new properties. For approximation on the canonical unit interval I := [-1,1], we will make our statements precise 
in the following subsection. 


3.1. One dimension 

Let K be the field of K or C. A function / : I —> /T is of bounded total variation if: 

Vf^ {\f{z)\dz<+^. (26) 

Ji 

Chebyshev polynomials of the first kind are defined by 1471 : 

T'„(x) = cos(ncos^^(.r)), for n e No, and .r e I. (27) 

A Chebyshev interpolant to a continuous function / : I —> is the approximation 

N-\ 

Pn{x) ^Y^c„T„(x), .xel, (28) 

n=0 


which interpolates / at the Chebyshev points of the first kind: 


PNix„) = f{x„) where 


cos 


2n + l 
2N 


for n - 0,... ,N - 1. 


(29) 


The Chebyshev basis has fast transforms between values at Chebyshev points and coefficients via fast imple¬ 
mentations of the discrete cosine transforms (DCTs). The (orthogonal) Chebyshev polynomials satisfy a three-term 
recurrence relation that can be used in Clenshaw’s algorithm m for 0{n) evaluation of interpolants. Compared with 
the best polynomial approximants, Chebyshev interpolants are near-best in the sense that their Lebesgue constants 
exhibit similar logarithmic growth. 
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Theorem 4 (Battles and Trefethen Let f be a continuous function on I, its N-point polynomial interpolant 

in the Chebyshev points of the first kind and p'^ its best degree-N — 1 polynomial approximation. Then: 

1 . 11/ - PnWo. < (2 + 2 log - l) 11/ - p*|U; 

2. if f has a derivative in I of bounded variation for some k > 1, \\f - Pa^||oo = 0(N^^) as N oo; and, 

3. if f is analytic in a neighbourhood of I, \\f — PatIIoo = 0(C^) as N —* oo for some C < 1; in particular we may 
take C — l/(M + m) if f is analytic in the closed Bernstein ellipse with foci +1 and semimajor and semiminor 
axis lengths M > I and m > 0. 

An interpolant can be constructed to any relative or absolute tolerance e by successively doubling the number of 
interpolation conditions, transforming values to coefficients, and determining an acceptable degre^ 

3.2. Two dimensions 

Numerous methods have been devised to approximate functions in more than one dimension. The straightforward 
generalization of the one-dimensional approach is to sample the function on a tensor of one-dimensional interpolation 
points and to adaptively truncate coefficients below a certain threshold. 

Consider the function f -.1^ ^ K, whose two-dimensional Chebyshev interpolant takes the form: 


m-1 n-l 

AijTfx)Tj(y). (30) 

i=0 j=0 

While the tensor approach in general suffers from the curse of dimensionality, it can still be competitive in two 
dimensions, scaling with 0(mn) function samples and 0(min(mn log n, nm log m)) arithmetic via fast two-dimensional 
transforms. 

The singular value decomposition of an m x n matrix A over K is the factorization ll53l : 

A = UEV*, (31) 

where U is an m x m unitary matrix over K, E is an m x n diagonal matrix of non-negative singular values, and V* 
is an n X n unitary matrix over K. The singular value decomposition reveals the rank of a matrix as the number of 
nonzero singular values. 

If we perform the singular value decomposition of the matrix of coefficients in ( [30l ), the approximation to / can 
be re-expressed as: 

k 

Psvoix, y) = ^ o-iUi{x)v\iy), (32) 

1=1 

where cr, are the singular values, and uf x) and v*(y) are univariate Chebyshev approximants with coefficients from 
the columns of U and the rows of V*, respectively, and where A is of rank k. It follows that psvo is the best rank-k 
approximant in L^(I^) to / that can be obtained for the original two-dimensional interpolant. For any given tolerance 
e > 0, a function / has numerical rank k^ if ll54l 

k, = inf|inf||/-/,|U<e||/|u|, (33) 

[ft J 

where the inner infimum is taken over all rank-k functions. 

Definition 5 (Townsend Definition 3.1 1|54|| ). For some e > 0, let be the numerical rank of f : ^ K, and 

and be the maximal degrees of the univariate approximations in the x and y variables. If k,;(m,; -H n^) < m^n^, we 
say the function / is numerically of low rank, and if k^ xi mm(m^, nf), then the function / is numerically of full rank. 


* This heuristic determination is usually based on, among other things, the relative and absolute magnitudes of initial and final coefficients, the 
decay rate of the coefficients, an estimate of the condition number of the function, and an estimate of the Lebesgue constant for a given degree. 
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A particularly attractive scheme for calculating low rank approximation in two dimensions can be described as 
a continuous analogue of Gaussian elimination ll54l and is a direct extension of the greedy algorithm in one dimen¬ 
sion ll48l Chapter 5]. This algorithm is studied in depth in Townsend’s DPhil thesis and implementations are found 
in Chebfun 1551 and ApproxFun. jl ES). In this algorithm, the function is initially sampled on a grid to locate its 
approximate absolute maximum. Two one-dimensional approximations are created in the x and y variables to inter¬ 
polate the function along the row and column that intersect at the approximate absolute maximum. After subtracting 
this rank-one approximation, the algorithm continues its search for the next approximate absolute maximum. After k 
iterations, it is clear that the approximant 


k 


Pov.(x,y) = '^Ai{x)Bi{y), 
?=1 


(34) 


coincides with / in the k rows and columns whose intersections coincide with an iteration’s approximate absolute 
maximum. As the size of the sampling grid increases, the approximate absolute maxima will converge to the true 
absolute maxima and in this sense we reproduce close aproximations to psvo- In terms of the degrees of the one¬ 
dimensional approximations m, n and the rank k, the algorithm scales with a search over 0{mn) function samples and 
0{k (m log m + n log ri)) arithmetic via fast one-dimensional transforms. 


Definition 6 (Townsend Definition 4.11 1|541 ). The function / ; —» AT is Hermitian if it satisfies the conjugate 

symmetry f(x,y) = f*(y, x) and it is non-negative definite, i.e.: 


for all a(x) e C(I). 


f£ 


a*(y) f(y, x)a{x) dy dx > 0, 


(35) 


When a bivariate function is Hermitian, even further savings can be obtained by drawing the analogy to the 
Cholesky factorization of a Hermitian matrix 15^ : 

k 

moiesky(x,y) = Y^Ai(x)A*(y). (36) 

1-1 

In this case, it is known that the function’s absolute maxima after every iteration are on the diagonal line y - x, leading 
to a reduction in the dimension of the search space. In addition, as they are conjugates only either the row or column 
slices may be computed and stored. 


3.3. An algorithm to extract the splitting of a fundamental solution 

Accurate numerical evaluation of a fundamental solution on or near the singular diagonal may not always be 
possible or may be more expensive ED. To avoid the numerical problems associated with the singular diagonal, we 
use Chebyshev points of the first kind in one direction and Chebyshev points of the second kind Il47l in the other 
direction. This ensures that the diagonal is never sampled. In terms of the DCTs, taking 2" points of the first kind is 
optimal and taking 2" -H 1 points of the second kind is nearly optimal. 

When both A(x, y) and B{x, y) in ( |25] l are not known a priori, but the fundamental solution itself can be evaluated, 
we can use such skewed grids in combination with the Riemann function to: 

1. approximate A(x, y) s -^iH(xi + ix 2 , Xi - ix 2 ,yi -i- iy 2 ,y\ - %); and subsequently, 

2. approximate the difference B(x, y) s <[)(x, y) - A(x, y) log |x - y|. 


4. The ultraspherical spectral method 

The ultraspherical spectral method of Olver and Townsend ifTSl represents solutions of linear ordinary differential 
equations of the form 


Jiu — /, IBu - c, 
1 


( 37 ) 


where ^ is a linear operator of the form 

d 

a^ix)-^ + ■ ■ ■ + a\_{x)—+aoix), (38) 

and S contains N linear functionals. Typically, S encodes boundary conditions such as Dirichlet or Neumann condi¬ 
tions. We consider u{x) in its Chebyshev expansion 


OO 

uix) = ^ u„Tn{x), (39) 

«-0 

so that u(x) can be identified by a vector of its Chebyshev coefficients u = (mq, mi, .. .)^. 

To solve such a problem efficiently, a change of basis occurs for each order of spectral differentiation, using the 
formula; 

d^T„{x)^( 0, Q<n<A-l, 

dx'^ \ 2^-\A-l)\nCl^!/x), n>A, ^ ^ 

where represents the ultraspherical polynomial of integral order A and of degree n. This sparse differentiation has 
the operator representation; 


/ A times 


Da = 2-^-‘(T- 1)! 


0 A 


A+l 


A + 2 


d>l, 


(41) 


and maps the Chebyshev coefficients to the d* order ultraspherical coefficients. 

Since in ( |38| l, each derivative maps to a different ultraspherical basis, the sparse differentiation operators are 
accompanied by sparse conversion operators such that ^ can be expressed completely in the basis of highest order N: 


So^ 


(1 0 

1 

2 




(1 0 

A 

A+2 


\ 

1 

0 

1 



A 

0 

A 


2 

2 



A+l 

A+3 



1 

2 

0 


Sa = 


A 

A+2 

0 


V 



‘ ‘) 


\ 



' •/ 


/!> 1 . 


(42) 


Here, >So maps the Chebyshev coefficients to the first order ultraspherical coefficients and maps the d* order 
ultraspherical coefficients to the (d-H 1)* order ultraspherical coefficients. Therefore, the conversion and differentiation 
operators can be combined in ^ as follows; 

(unDn + H- ■ ■ ■ -H goSn-i ■ ■ ■ So) u = *5^-1 • ■ ■ *Sof, (43) 


where u and f are vectors of Chebyshev expansion coefficients. Were the coefficients ai(x), i - 0,...,N, constant, 
then ( |43| would represent a linear recurrence relation in the coefficients u of length at most 2N + 1. However, the 
coefficients are in general not constants, so the multiplication operators in Chebyshev and ultraspherical bases are also 
investigated in lfT3l . Let 

OO 

a(x) ^Y^anT„(x). (44) 

n=0 

Then it is shown in Ha that multiplication can be represented as a Toeplitz-plus-Hankel-plus-rank-one operator; 


Mo[a] = 


1 

2 


'2ao 

a\ 

02 

• • ■'1 


'0 

0 

0 

■ • 

ai 

2ao 

a\ 


+ 

fli 

02 

02 


a2 


2ao 



02 

0.2 

O^ 


\ 



‘ •} 


\ • 



• ‘ / 


(45) 
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For /! > 0, an explicit formula for the entries is given in El and a three-term recurrence relation is shown in 1541 
Chap. 6]. By the associative and distributive properties of multiplication, the recurrence relation for the multiplication 
operators is derived from the recurrence relation for the ultraspherical polynomials: 




n + 1 


(46) 


Since we assume the coefficients a,(x) to be continuous functions with bounded variation on I, let m denote the highest 
degree Chebyshev expansion such that for some e > 0: 


m-I 

Oiix) - UinTnix) 
n=0 


< e||a,(.x)||oo 


for i -0,...,N. 


(47) 


Then in this way, the system 


Su = c, 

{MN[aN]I>N + MN[aN-\]SN-\I>N-\ H-+ AlA'[ao]‘Siv-i ■ ■ ■>So)u = Sn-\ ■ ■ ■>Sof, (48) 

is almost banded with bandwidth 0{m). The proposed OiirP'n) solution process for such systems is the adaptive QR 
factorization, generalizing (F. W. J.) Giver’s algorithm for second-order difference equations l58ll . In this factorization, 
the forward error is estimated at every step in the infinite-dimensional upper-triangularization to adaptively determine 
the minimal order n required to resolve the solution below a pre-determined accuracy. Since the unitary transforma¬ 
tions implied by Q preserve the rank structure, the back substitution is also performed with 0{m^n) complexity. 

Figure [T] shows the typical structure of the system and an example of the type of singularly perturbed boundary 
value problem that it can solve efficiently. 




Figure 1: Solution of e(e + x^)u''(x) = x u(x), u(-l) = 1, h( 1) = 0 via the ultraspherical spectral method. Left: the structure of the system. Right: a 
plot of the solution for e = 10^“*. In this case, a Chebyshev expansion of degree 3,276 is required to approximate the solution to double precision. 


4.7. Almost-banded spectral methods in other bases 

The key elements of the ultraspherical spectral method are a graded set of bases that permit banded differentiation 
and conversion within the set of bases, and multiplication operators for variable coefficients. Other examples where 
a graded basis can be exploited are the Jacobi polynomials (which include Legendre and ultraspherical polynomials 
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as special cases), and the generalized Laguerre polynomials. Hermite polynomials, which form an Appell sequence, 
satisfy - 2n//„_i(x), and therefore do not require other bases for conversion. 

From the three-term recurrence relation satisfied by orthogonal polynomials ||59|: 


X7r„(x) = a„7r„+,(x) + j5„nn(x) + y„7r„_i(x), 
it follows that multiplication by x is tridiagonal; 


(49) 


M[x] = 


71 y6l Qfl 


72 Pi Q'2 


(50) 


Therefore, banded multiplication operators in orthogonal bases can be derived from the recurrence relation: 

n > 1. 


\ j 


(51) 


and consequently variable coefficients represented as interpolants have a finite-bandwidth operator form. To numeri¬ 
cally determine such variable coefficients practically requires fast transforms. Among the many possibilities, see 
for a new approach for a fast FFT-based discrete Legendre transform. 


5. Ultraspherical spectral method for singular integral equations 

In the following definitions, we identify C with and let F be bounded in C. 

Definition 7 (Kress IIH). A real- or complex-, scalar- or vector-valued function / defined on F is called uniformly 
Holder continuous with Holder exponent 0 < a < 1 if there exists a constant C such that 

l/(x)-/(y)|<C|x-yr, for x,yeF. (52) 

By C’'’“'(F) we denote the space of all bounded and uniformly Holder continuous functions with exponent a. For 
vectors, we take | ■ | to be the Euclidean distance. With the norm 

iirii ir/ M , l/(x)-/(y)l 

||/||o,a. sup |/(x)| H- sup —^--—, (53) 

xer X.yer |X - y| 

XTty 


the Holder space is a Banach space, and we can further introduce C^’"(F) as the space of all differentiable functions 
whose gradient belongs to C’'’“'(F). 

Definition 8 . Let / e C°’“(F). The Cauchy transform over F is defined as: 


Cr/(z):=^ r for z e C \ F. 

27n Jr ^-z 

The Cauchy transform can be extended to z e F with integration understood as the Cauchy principal value. 
Definition 9. Let / e C®’“(F). The Hilbert transform over F is defined as: 


(54) 


'Hr/(z):=-f^d^, for z e F, 
where the integral is understood as the Cauchy principal value: 

/^d,= ilimr 

C-z ' np-^o Jr\rte;p) (-z 



(55) 


( 56 ) 


where F(z;p) 6 F : |^ - z| < p). 
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Lemma 10 (Sokhotski-Plemelj EH ED). Iff e C“’“(r), then: 


■Hr/(z) = i[C^+C-]/(z), 


( 57 ) 


where denotes the limit from the left/right ofY. 

With the Hilbert and Cauchy transforms, further integrals with singularities can be defined. 
Definition 11. For / e C°(r) the log transform over F is defined as: 


Xr/Cz) - flog 1^ - zm) for z € C. (58) 

?!■ Jr 

For / 6 C''“(F) the derivative of the Hilbert transform is defined as: 

'Wf/(z):=-f for ^ ^ T- (59) 

TrJr - zY 

where the integral is understood as the Hadamard finite-part ll^ 1^ : 


If no 

^Jr - zY 




- lim 

TT p-»0 


r f(o 

Jr\r(z;p) (^ “ zY 


d^- 



where F(z; p) is defined as in Definition]^ 


Remarks. 


(60) 


1. The Sokhotski-Plemelj lemma offers a convenient way to compute the Hilbert transform via the limit of two 
Cauchy transforms. 

2. The use of the Cauchy principal value and the Hadamard finite-part allows for the regularization of singular and 
hypersingular integral operators, respectively. 

On a contour F, we expand the kernel of the singular integral equation (0 in the following way: 


iHu - f, Su — c, 


( 61 ) 


for 

-I- log |y - x\KYx, y) -H K^x, y)] M(y) dy, 

where Ki, K 2 , and K 4 are known continuous bivariate kernels, / is continuous, !B contains N linear functionals, 
and u is the unknown solution. If in ( |6T] l, we replace the bivariate kernels with low rank approximations, 

kx 

Kx{x,y) :=i'^AxYx)B 2 ,i(y), for T=l,2,3,4, (62) 

1=1 

we achieve at once two remarkable things: firstly, the approximations are compressed representations of the kernels; 
and secondly, the separation of variables in the low rank approximation allows for the singular integral operators to 
be constructed via the Definitions 1^ and [TTl 

In the following two subsections, we consider the case where F is the unit interval, and emulate the construction 
of the ultraspherical spectral method for ODEs to arrive at an almost-banded system to represent ( |6T] i. In this setting, 
we must use weighted Chebyshev bases to accomplish this task. Note that alternative spectral methods for open arcs 
are discussed in ED. 
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5.1. Inverse square root endpoint singularities 

Indeed, the Hilbert transform of weighted Chebyshev polynomials is known 




(- 1 . 1 ) 


Tnix) 

Vl - 


0 , n = 0, 
Cll\(x), n>l, 


( 63 ) 


This operation can then be expressed as the banded operator from the weighted Chebyshev coefficients to the ultras- 
pherical coefficients of order 1; 

(0 1 

1 

1 . (64) 


n 


(-1,1) - 

Upon integration with respect to x, we obtain an expression for the log transform: 

T„(x) 


-C(-i 


(- 1 , 1 ) 


Vl — 


- log 2, n - 0, 
Tn(x) 

-, n > 1, 


(65) 


or as an operator from the weighted Chebyshev coefficients to the Chebyshev coefficients: 

-log 2 

^(-U) = -i 


( 66 ) 


In addition, upon differentiation with respect to x, we also obtain an expression for the derivative of the Hilbert 
transform: 


^;i,i) 


Tn(x) 


Vl - X^ 


0 , 


n = 0,1, 


C%(x), n>2, 


(67) 


This operation can then be expressed as the banded operator from the weighted Chebyshev coefficients to the ultras- 
pherical coefficients of order 2: 

p 0 1 

1 

- 1 ■ ( 68 ) 

. ‘ ■/ 

Lastly, the orthogonality of the Chebyshev polynomials immediately yields for the functional 

2r/:= - r fiOdi 

^ Jr 

T„(x) 


(69) 


the following: 


K- 1 , 1 ) 


Vl - x^ 

or as a compact functional on the weighted Chebyshev coefficients: 


1, n = 0, 

0 , n > 1, 


5:(_u) = (i 0 0 ■■■). 


(70) 


(71) 
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Combining the integral operators together with the bivariate approximations, we define; 


h 

:= X (72) 

(=1 

^=2 

^(-x.x)[K2\ := 2 MdA2Ax)YHi.-U)MQ[B2,i(y)l (73) 

(=1 
^=3 

:= 2 Ado[A3,,(^)]X(-i,i)A1o[B3,,(y)], (74) 

k4 

S(-i,i)[-^^4] := ^ Ado[A4,i(-r)]2(_ip)A1o[7J4,iCv)]- (75) 

!=1 

Then, we can reduce singular integral equations of the form ( |6T] l into an infinite-dimensional almost-banded system: 

Su = c, 

+ S\H{-\^\)[K2\ -I-.SnSo(-£(-i,i)[fif3] -H2(_i4)[^r4]))u = *SnSof. (76) 


This system can be solved directly using the framework of infinite-dimensional linear algebra lfT4l . built out of the 
adaptive QR factorization introduced in II3. 


5.2. Square root endpoint singularities 

The Hilbert transform of weighted Chebyshev polynomials of the second kind is also known 


'Hi [Un(x)^l - x^\ = -7’„+i(x), n>0. 


(77) 


This operation can then be expressed as the banded operator from the weighted ultraspherical coefficients of order 1 
to the Chebyshev coefficients: 

0 

-1 

-1 . (78) 


Upon integration with respect to x, we obtain an expression for the log transform: 


£l [Unix) Vl - X^] = 


-^log2-H i7’2(;c), n = 0. 


1 / T„+2ix) Tnix) 


(79) 


, n > 1, 


2 \ n -I- 2 i 

or as an operator from the weighted ultraspherical coefficients of order 1 to the Chebyshev coefficients: 


-Ci = 


^-5 log 2 

0 

1 

4 


_\ 

2 

0 4 


(80) 


In addition, upon differentiation with respect to x, we also obtain an expression for the derivative of the Hilbert 
transform; 

[Unix) Vl - x2j = -in + 1)C''V-^), « > 0, (81) 
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This operation can then be expressed as the banded operator from the weighted ultraspherical coefficients of order 1 
to the ultraspherical coefficients of order 1 : 


(-1 




-2 


-3 


(82) 


Lastly, the orthogonality of the Chebyshev polynomials of the second kind immediately yields for Sj: 

Si[t/„(^)Vr^] = | 2’ ” = (83) 

or as a compact functional on the weighted Chebyshev basis: 

= 0 0 ■■■). (84) 

Combining the integral operators together with the bivariate approximations, we define: 


h 


:= J]Mi[Ai,,(^)]'?fi'Afi[Bi,,( 3 ;)], 

(85) 

kl 

mKi] :=J]MQ[A2jix)]9{iMi[B2j(y)], 

i=l 

( 86 ) 

h 

:= ^Mo[A3,,(^)]i:iAfi[B3,i(y)], 

/=! 

(87) 

h 

( 88 ) 


!=I 


and in the framework of infinite-dimensional linear algebra IIT4ll . we may solve singular integral equations of the 
form (|6T]) via the almost-banded system: 


Su = c, 

{9{{[Ki] + SoiniKi] + 2:i[/:3] + MK4]))u = *Sof. 


(89) 


Let nix + 'Wy denote the largest sum of degrees of the bivariate Chebyshev expansions of the integral kernels such 
that for some e > 0 : 


KAix,y) - Y^AA,iix)BAj(y) 


i=l 


< e\\KA(x,y)\\o 


for 4 = 1 , 2 ,3,4. 


(90) 


Then, the complexity of the adaptive QR factorization is 0{{mx + my)^n) operations, where n is degree of the resulting 
weighted Chebyshev expansion of the solution. This is reduced to Oiinix -i- my)n) operations by pre-caching the QR 
factorization. 


Remarks. 

1. The observ ation that |d^| = d^ on I allows us to relate fine integral formulations with the operators of Defini¬ 
tions!^ and [00 

2. Mixed equations involving derivatives and singular integral operators are also covered in this framework. 

3. It is straightforward to obtain the singular integral operators on arbitrary (complex) intervals {a,b) using an 
affine map. 


^Both variants of the singular integral operators are implemented in SingularIntegralEquations. j 1. 
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5.3. Multiple disjoint contours 

Singular integral equations on a union of disjoint intervals T = Ti U r 2 U ■ ■ • U are covered in this framework. 
We can decompose (|M]l as 


f Si 
-^1,1 
-^2.1 

S2 ■ 
-^1.2 ■ 

-^2,2 ■ 

■ Sd] 

-^l.d 

^2,d 

'Ui' 

U2 

— 

'c' 

fl 

f2 

(91) 


^d.2 ■ 

■ ^d,dj 

Md. 





where each S, is a set of linear functionals and Jlij - Jlrilrj- The diagonal blocks are equivalent to the previous 
case considered, hence result in banded representations. The off-diagonal blocks can be constructed directly by 
expanding the entire non-singular kernel in low rank form and using the compact functionals 2(_i i) or Ei. The 
resulting representation is, in fact, finite-dimensional and hence every block is banded. 

Here, we show how a block-almost-banded infinite-dimensional system can be interlaced to be re-written as a 
single infinite-dimensional and almost-banded system. Re-ordering both vectors (ui, U 2 ,..., and (fi, f 2 ,..., f^/)^ 
to: 


U = (uifl M2,0 ■ ■ ■ Ud ,0 Ml,l “2,1 ■ ■ ■ “d.l ■ ■ ■) (92) 

F = (/i,o h.o ■ ■ ■ fd,o fi,i fi.i ■■■ fd,\ ■■ ■) (93) 

amounts to a permutation of almost every row and column in (|9T|l. Define each entry of and 21 by: 


95,',2 - S|(/_l)niodrf)+l,L^J,p 


(94) 

(95) 


where the last two indices in each term on the right-hand sides denote the entries of the functional or operator. This 
perfect shuffle allows for the system (|9T]l to be re-written as the almost-banded system 



(96) 


5.4. Diagonal preconditioners for compactness 

We now show that our formulations leads to equations whose operators are compact perturbations of the identity. 
For well-posed (integral) equations, this ensures convergence Ha. We show this for the singular operators in equa¬ 
tions ([TT} and ( [T^ defined on the canonical unit interval and in suitably chosen spaces. Note that a similar analysis is 
performed in ||2TI . Since we are working in coefficient space, we consider the problem as defined in spaces. In the 
case of Chebyshev expansions, this corresponds to Sobolev spaces of the transformed function m(cos ff). 

Definition 12 (Olver and Townsend llT3l ). The space c C“ is defined as the Banach space with norm: 


Hull,. = 


00 


A 


^|M,|2(k+ 1)2-1 

k=0 


< 00. 


(97) 


Let P„ = (4,0) be the projection operator. 


Lemma 13. For the Dirichlet problem singular integral operator in o, if O takes the form ( |25| l with A and B 
analytic in both x and y and if we take 7? to be 


( 1 

log 2 


■R = 2 


1 


2 


3 





(98) 


V 
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then 

{£(-i,i)[7tA] + 2(_i,i)[;rB])?? = / + 7C, (99) 

where 'K : is compact for /I 6 M. 

Proof. Since A(x, x) = -(2ny^, we let A(x,y) = A(x,y) - A(x, x) and separate the operator 0 as: 

X(_i,i)[7rA(x,x)] +X(_i,i)[7rA(x,3;)] + S(_i,i)[7r5]. (100) 

It is straightforward to show 

£^^iy^[jtA{x,xm^I:£l^{i ( 101 ) 

Then, we need to show that the remainder is compact. Since; 


- -C(-i,i)|| ^0 as n ^ oo, 


( 102 ) 


-C(-i,i) \ t\ ^ t\ is compact. Compactness of is implied by its finite-rank. Expanding A and B in low rank 

Chebyshev approximants, we have: 


f ki 


kp 


^ Ado[Ai,,(^)]i:(-i,i)Mo[A2,,(y)] + ^ Mo[BUxm-i,i)Mo[B2j(y)] 


\ i=l i=l 

Since A and B are analytic with respect to y, then for every i and for every A 6 ! 


'R. 


M)[A 2 ,,(y)] : 
Mo[B2,m ■ tl-x 


MQ[A2,i(ym : (I ^ (I 
Mo[B 2 ,i(y)YJ< ■■ ^ 


(103) 


(104a) 

(104b) 


are bounded. Compactness follows from the linear combination of a product of bounded and compact operators being 
compact. □ 

Lemma 14. For the Neumann problem singular integral operator in 0, if O takes the form ( |25| l with A and B 
analytic in both x and y and if we take R to be 


(1 


1 

2 


R ^-2 


V 





• 




■) 


(105) 


then: 

where the two primes indicate: 


(■//'[-ttA] + Xi[7rA"] + l.i[7tB''])R ^I + 'K, 


A"ix,y) 


d'^Aix, y) 


dx 2 dy 2 


x,y=(x,0),(y,0) 


and where “TC : ^ is compact for 4 € K. 

Proof Since A(x, x) = -(27r)“', we let A{x,y) = A{x,y) - A( 2 i:, x) and separate the operator 0 as: 

'H{[-7tA{x, x)] + 'H;[-7TA(x,y)] + Li[nA"] + Y-AnB"]. 


It is straightforward to show: 

'H[[-nA{x, x)]R - I: 


(106) 

(107) 


(108) 

(109) 
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Then, we need to show that the remainder is compact. Since: 


Wn'R'Pl - ‘??ll ^ 0 as n ^ oo, (110) 

H ^ compact. Furthermore, showing boundedness of .So,-Ci[,Si : ^ and 'TTj : is 

straightforward. Expanding A, A" and B" in low rank Chebyshev and ultraspherical approximants, we have: 


TT - 2 Ml [A,,i{x)YH{Mx[A2j(y)] 

. i=l 

+ *So 2 Mo[A",(x)]XiMi [A'lfy)] + 2 Mo[B'i',,(x)]SiMi [B'l,{y)] 

,;=1 i=l 


n. 


Since A and B are analytic with respect to y, then for every i and for every d 6 M: 


MAAijiy)] : el 


^A+\ 


H'Mi[A2,,0')] : el ^ el. 


( 111 ) 


( 112 ) 


are bounded. Compactness follows from the linear combination of a product of bounded and compact operators being 
compact. □ 

Remarks. 


1. For complicated fundamental solutions whose bivariate low rank Chebyshev approximants have large degrees, 
preconditioners such as those in Lemmas 13 and 14 allow for continuous Krylov subspace methods or conjugate 
gradients on the normal equations to converge in a relatively fewer number of iterations compared with the un¬ 
preconditioned operators. Furthermore, the low rank Chebyshev approximants allow for the operator-function 
product to be carried out in 0 {{m + n) log(m -H n)), where m is the largest degree of a multiplication operator 
and n is the degree of the Chebyshev approximant of the solution. Iterative solvers are outside the scope of this 
article, however. 

2. Operator preconditioners can also be derived which would yield similar I + % results. However, working 
in coefficient space allows for a simpler exposition. 


5.5. Numerical evaluation of Cauchy and log transforms on intervals 

Fast and spectrally accurate numerical evaluation of the scattered far-held can be derived from Clenshaw-Curtis 
integration of the fundamental solution multiplied by the density. For each evaluation point, the fundamental solution 
can be sampled at the 2N roots of the 2A^* degree Chebyshev polynomial, where N is the length of the polynomial 
representation of the density. Since the resulting density may be as complicatecQ as the fundamental solution itself, 
doubling the length is sufhcient to resolve the coefficients of the fundamental solution multiplied by the density. 

It is well known that such an evaluation technique is inaccurate near the boundary lICTl . In the context of Riemann- 
Hilbert problems, spectrally accurate evaluation near and far from F can be obtained by exact integration of a modihed 
Chebyshev series that encodes vanishing conditions at the endpoints. 

Consider the modified Chebyshev series: 

fo(x):=l, fi(x):=x, t„{x) T„(x) - T„_ 2 (x), n>2. (113) 

If we expand m in a Chebyshev series and this modified Chebyshev series: 

CO OO 

u(x) = ^ UnT„{x) = ^ U,J„{x), (114) 

n=0 n=0 


^In the Helmholtz equation, for example, both the density and the fundamental solution are oscillatory with the same wavenumber. 
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then we have the relation: 


'uf 


o 

1 


'mo' 

U\ 


1 0 -1 


Ml 

U2 

— 

1 0 -1 


U2 

s • . 




K • Z 


(115) 


Therefore, any finite sequences and can be transformed to the other in 0{N) operations, either via 

forward application of the banded operator, or via an in-place back substitution. 

Lemmas 16 and 17 contain formulae for Cauchy transforms of weighted Chebyshev polynomials evaluated in 
the complex plane. These were originally derived in this form in Il68l . based on results in ll69l iTOl 1711 1721 . These 
formulae are adapted in Lemma 18 for the log transform as well. 


Definition 15. Define the Joukowsky transform: 

(116) 

and one of its inverses: 

:=z- V^V^, (117) 

which maps the slit plane C \ I to the unit disk. 


The Joukowsky transform is useful for proving and summarizing the following results. 

Lemma 16 (Lemma 5.6 1^ ). For k > 0.' 

Ci[Vl-o2t/,](z) = (118) 

Proof. We verify that the Sokhotski-Plemelj lemma is satisfied. Note that for x - cos 6 we have: 

lim77^(.r + ie) - x + i V1 - = cos0 + isin0 = e’"’®. (119) 

e\0 


It follows that: 


7, ^(cos 6 + - 7,'(cos 9 - e - Q<k+i)e 

lim--> 

f\o 2i 2i 


= - sin(k -H 1)0 = -Ukicos 6) Vl - cos^ 9. 


Lemma 17 (Lemma 5.11 Il68l l. For k > 2: 

Q- 1 , 1 ) 

C(_i,i) 

Q”i,i) 


1 


Vl — <>2 

o 


Vl - <>2 
fk 


( 2 ) = 


(z) = 


2 Vz - 1 Vz + 1 


IZ 


2 Vz - 1VF-TT 2 
(z) = -i7;'(z)'^-'. 


- -, and 


Vl - <>2 

Proof. The first two parts follow immediately from the Sokhotski-Plemelj lemma. The last part follows since: 

cos k9 - cos(k - 2)9 


sin(A: - 1)0 = 


2 sin0 


( 120 ) 

□ 

( 121 ) 

( 122 ) 

(123) 

(124) 

□ 
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We extend these results here to the log transform. 


Lemma 18. 




Xi[t/,Vl-o2](z)= -!R 
1 


n\zf 


k + 2 


■£(- 1 . 1 ) 

-£(- 1 , 1 ) 

-£(- 1 , 1 ) 

-£(- 1 , 1 ) 


. Vi -o2 

o 


Vl 

fl 


Vl — 
fk 


(z) = -log|7+‘(z)| -log2, 
(z) = -%Jl\z), 

(z) = log|7;'a)|+log2-5l 


n\zf 


and 


Vl 




n\zf~^ n\zf 


k-2 


(125) 

(126) 

(127) 

(128) 

(129) 

(130) 


Proof. These formulae follow from integrating the formulae for Cauchy transforms and taking the real part. We can 
compute the indefinite integrals directly lIMl §5.4.4]; 


r . - ^ I — = dz = - log 7+'(z), 

J Vz^VzTT 

1 -z . .-I, 


f 


Vz- 1 Vz -t-1 


dz = 7;‘(z), 


2/4;. 

:/ 4 ;. 




(z) dz = 


log 7+'(z), and 


^k^._ 7;‘(z)'-* 


(z)'=dz = 


k+1 


k- 1 


for k>2. 


We also have the normalization for z ^ -t-oo: 


J f(x) log(z - x) Ax ^ log z ^ f(x)Ax + 0{z^). 


Note that: 


hence; 


1 


7+ (z) ~ ;:r + £^(^ ^ ^ 

2z 

log74'(z) =-logz-log2-i-(9(z^') as z- 


(131) 

(132) 

(133) 

(134) 

(135) 

(136) 

(137) 

□ 


These formulae can be generalized to other intervals, including in the complex plane, by using a straightforward 
change of variables: 


r f, r 

Ma.b)J(Z) - ----L(-l,l) 


/ 


b + a b — a 

-+-o 

2 2 

\b — a\ _^_^\b — a\ 

2 n 


b + a - 2z 


log 


i:a 


b — a 
b + a b — a 


x\ Ax. 


(138) 
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6. Applications 

6.1. The Faraday cage 

The Faraday cage effect describes how a wire mesh can reduce the strength of the electric field within its con¬ 
finement. This phenomenon was described as early as 1755 by Franklin ITT^ §2-18] and in 1836 by Faraday 11741 . 
While the description of the phenomenon is quite prevalent in undergraduate material on electrostatics, a standard 
mathematical analysis has been missing until only recently by Martin ITTSI and Chapman, Hewett and Trefethen llT^ . 
In IT^ . three different approaches are considered for numerical simulations: a collocated least squares direct numeri¬ 
cal calculation, a homogenized approximation via coupling of the solutions at multiple scales, and an approximation 
by point charges determined by minimizing a quadratic energy functional. 

In IIT^ . it is shown that the shielding of a Faraday cage of circular wires centred at the roots of unity is a linear phe¬ 
nomenon instead of providing exponential shielding as the number of wires tends to infinity for geometrically feasible 
radii, i.e. radii that prevent overlapping. In their synopsis, it is claimed that a Faraday cage with any arbitrarily shaped 
objects will not provide considerably different shielding in the asymptotic limit. Here, we confirm this observation 
with infinitesimally thin plates of the same electrostatic capacity as wire^angled normal to the vector from the origin 
to their centres. Our numerical results are in excellent asymptotic agreement with those presented in iTh). Departing 
from the practical case of normal plates, we also consider infinitesimally thin plates angled tangential to the vector 
from the origin to their centres. In this case, we escape the practical material limit on the number of shields as an 


infinite number of plates can be modelled independent of radial parameter. 

We seek to find the solution to the Laplace equation such that, in addition: 

Am(x) = 0, for X € O, (139a) 

m(x) = Mo, for X e F, (139b) 

m(x) = log |x - y| -I- (9(1), as |x - y| —> 0, (139c) 

m(x) = log |x| -H o(l), as |x| ^ oo. (139d) 

Since this is a Dirichlet problem, we begin by splitting the solution u - u' + u\ where: 

m'(x) = log |x - y| = 2;ril)(x, y), (140) 


is the source term with strength In located at y = (2,0), as in Ea. We represent in terms of a density with the 
single-layer potential equal to the effect of the logarithmic source. Alone, this represents a solution to the Laplace 
equation with Dirichlet boundary conditions on F. To satisfy condition ( |139b] i, we augment our system to ensure there 
is a constant charge of zero on the wires and plates: 

dF(y) = 0, (141) 

though each wire may individually carry a different charge, and the unknown constant mq to accommodate this condi¬ 
tion. Figure]^ shows the numerical results for shielding by normal and tangential plates. Figure]^ shows a plot of the 
convergence of the density coefficients and the field strength at the origin for various parameter values. 



6.2. Helmholtz equation with Neumann boundary conditions 

The mathematical treatment of the scattering of time-harmonic acoustic waves by infinitely long sound-hard obsta¬ 
cles in three dimensions with simply-connected bounded cross-sections leads to the exterior problem for the Helmholtz 


equation: 

(A + k^)u{x) - 0, 

for k 6 K, X e Q, 

(142a) 


du(x) _ ^ 
dn(x) 

for X e F, 

(142b) 


ita = 

r^+oo y dr I 

for r |x|. 

(142c) 


'^This corresponds to plates of width 4r where r is the wire radius. 
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Figure 2: Left: a plot of the solution w(x) with 10 normal plates with radial parameter r = 10 ^. Right: a plot of the solution w(x) with 40 tangential 
plates with the same radial parameter, surpassing the material limit in the original numerical experiments ESI In both contour plots, 31 contours 
ai‘e linearly spaced between -2 and +1. 




Figure 3: Left: a plot of the Cauchy error of successive approximants for the solution of Laplace’s equation with 10 normal plates with r = 
10“^, corresponding to the left plot in Figure]^ where n is the total number of degrees of freedom. The + indicates where the adaptive QR 
factorization terminates in double precision, and with this approximation the forward error is ||wo + Sr[duydn] dr(y) - = 2.90 x 10“^^ and 

^r(y )||2 = 1-47 X 10“^^. Right: a plot of the field strength in the center of the cage versus the number of plates. The dashed lines 
represent results for normal plates, while the solid lines represent results for tangential plates of the same electrostatic capacity. The normal and 
tangential plates exhibit different asymptotic scalings. 
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Equation ( |142b| i enforces sound-hard obstacles, while equation ( |142c| i is the Sommerfeld radiation condition ifTTl . an 
explicit radiation condition at infinity. Consider an incident wave with wavenumber k and unit direction d: 

m'(x) = (143) 

We wish to find the scattered field m* such that the sum u - u' + u’‘ satisfies the Helmholtz equation in the exterior. 

The fundamental solution of the Helmholtz equation is proportional to the cylindrical Hankel function of the first 
kind of order zero IfTSl §8.405]: 

0(x,y)= (144) 

and the Riemann function is also well known BTII for the Helmholtz equation: 

m{z, ZO, ^o) = Mk ^J(z-zo)(^-m (145) 

Figure shows the rank structure of the bivariate kernels and the total solution with a set of randomly generated 
screens between [-3,3]. N.B. it is known that ll79ll collinear screens have reduced off-diagonal numerical ranks 
comparedwith randomly oriented screens. 



Figure 4: Acoustic scattering with Neumann boundary conditions from an incident wave with i = 100 and d = (1/ V2, -1/ V2). Left: a plot of the 
numerical ranks of Jo{k\x - y|) connecting domain i to domain j, where it can be seen that interaction between domains is relatively weaker than 
self-interaction. Right: a plot of the total solution. 1,392 degrees of freedom are required to represent the piecewise density in double precision. 


6.3. Gravity Helmholtz equation with Dirichlet boundary conditions 
The Helmholtz equation in a linearly stratified medium: 





(A -1- £ -1- X2)u(x) - 0, 

for £ e M, x 6 Q, 

(146a) 






m(x) = 0, 

for X € E, 

(146b) 
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du 
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dx2 = 0, 


(146e) 


models quantum particles of fixed energy in a uniform gravitational field ifSTl . Equation ( |146b| ) enforces sound-soft 
obstacles, while equations (|146c|i-(ri46eli form an explicit radiation condition at infinity derived in 1571 . 
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The fundamental solution of the Helmholtz equation in a linearly stratified medium is derived in lISOll : 




1 

Ttt 





12^ 


df 

t ' 


(147) 


Numerical evaluation via the trapezoidal rule lISTl along a contour of approximate steepest descent on the order of 10^ 
evaluations per second is reported in Ell. This equation is also known as the gravity Helmholtz equation. 

Consider an incident fundamental solution with energy E and source y: 


m'(x) = <l)(x,y). 


(148) 


We wish to find the scattered field m* such that the sum u - u' + satisfies the gravity Helmholtz equation in the 
exterior. In addition to the fundamental solution, we require the Riemann function of the PDO. With the prospect of 
deriving a fast numerical evaluation in future work, we prove the following theorem in [Appendix A| 

Theorem 19. The Riemann function of the gravity Helmholtz equation, where c(xi , X 2 ) — E + X 2 and therefore 
C(z,f) - j + ^ has the power series: 


in(z, zo, 4'o) = 1 + Xi z ~ 


(149) 


/=1 ;=1 

where the coefficients Ajj satisfy ( |A.3| )-( [A3] t, and the integral representation: 

I ny+ioo ^ 

2m Jy. 


, _I exp ISiE (( 5 ^ - 

y-ioo ys^ — ujAl 

y + (v - t/)5| 


(150) 


where 93(z, zo, ^o) = y{z - - ^o) where E - — 

4 81 

Figure 1^ shows the total solution to the gravity Helmholtz equation with Dirichlet boundary conditions and the 
2-norm condition number of the truncated and preconditioned system. 
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Figure 5: Acoustic scattering with Dirichlet boundary conditions from an incident fundamental solution (I>(x, y) with y = (0, -5) and £” = 20 against 
the sound-soft intervals ((-10, -3), (-5,0)) U ((-2,5), (2,5)) U ((5,0), (10, -3)). Left: a plot of the 2-norm condition number of the truncated and 
preconditioned system with n degrees of freedom. Right: a plot of the total solution. 332 degrees of freedom are required to represent the piecewise 
density in double precision. 
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6.4. Helmholtz equation with nearly singular Dirichlet boundary data 

In this application, we consider the Helmholtz equation with nearly singular Dirichlet boundary data. Consider the 
scattering of a collection of point sources arbitrarily close to a sound-soft obstacle. If we parameterize the locations 
of the point sources by the family of Bernstein ellipses Ep, then we know that the Chebyshev series representation 
of the incident wave will have degree which scales as « = (9((logp) ') as p —> 1. Additionally, as the point sources 
approach the boundary, the integral operator has bandwidth 0(k), independent of the Bernstein ellipse parameter. This 
is a challenging scenario for conventional integral equation solvers since a piecewise polynomial approximation to the 
nearly singular boundary data may not be much more efficient than a global representation. Furthermore, if the point 
sources are allowed to move freely on the Bernstein ellipse, then no adaptivity may be used to uniformly accelerate 
the solvers. The demonstrations in this section are also applicable to the important problem of many micro swimmers 
in Stokes flow approaching an obstacle, as the swimmers can be modelled as point sources, see ll82l . 

This set of problems completely demonstrates the scaling OHm^ + m.yf'n) of our algorithm: the bandwidth scales 
with the wavenumber, and the degree scales with the reciprocal of the log of the Bernstein ellipse parameter. Addi¬ 
tionally, a partial QR factorization of the singular integral operator may be cached or precomputecH resulting in the 
reduced 0{{mx + ;«,,)«) complexity for additional solves. Figure shows the scalings of the computation for three 
wavenumbers and varying Bernstein ellipse parameters. The figure also shows the solution of the Helmholtz equation 
with 100 nearby source terms. 




Figure 6: Timings to solve the Helmholtz equation (left) with nearly singular boundary data, and a sample of the solution (right). Left: timings ai‘e 
illustrated for three wavenumbers and multiple Bernstein ellipse parameters resulting in Chebyshev expansions reaching degrees on the order of 
half a million. The high data set show the partial QR factorization and the back substitution, while the low data set show the reduced time with a 
cached QR factorization. Right: a sample solution atk = 50, with 100 point sources uniformly distributed in angle on the top half of the Bernstein 
ellipse £'i.o 5 with charges +1 on the right half and -1 on the left half, resulting in the symmetric output. 


7. Numerical Discussion & Outlook 

The software package SingularIntegralEquations. jl IfTSl written in the Julia programming language ifThl 
El implements the banded singular integral operators, methods relating to bivariate function approximation and 
construction with diagonal singularities, fast & spectrally accurate numerical evaluation of scattered fields and several 
examples including those described in this work. Built on top of ApproxFun. j 1, SingularIntegralEquations. j 1 
uses the adaptive QR factorization described in lfT3l and acts as an extension to the framework for infinite-dimensional 
linear algebra. All numerical simulations are performed on a MacBook Pro with a 2.8 GHz Intel Core i7-4980HQ 


^The cached QR factorization can be adaptively grown without re-computing from scratch by exploiting the fact that the operator is banded 
below, thus the number of degrees of freedom (n) needed to resolve the solution within a prescribed tolerance need not be known apriori. This 
automatic caching of the QR factorization is implemented in ApproxFun. j 1. 
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processor and 16 GB of RAM. While timings are continuously being improved, Table [T] shows the current timings 
to solve the problems in section All the numerical problems relating to our applications have been abstracted so 
that to explore a new elliptic PDE in SingularIntegralEquations. j 1, the user only needs a fast evaluation of the 
fundamental solution and its Riemann function. 


Table 1: Calculation times in seconds to solve the problems in section]^ Evaluation of the scattered field is reported per target. Timings for the 
Laplace equation are for 10 normal plates. 



Kernel assembly 

Adaptive QR 

Evaluation of scattered field 

Faplace 

0.888 

0.518 

0.0000135 

Helmholtz {k = 100) 

1.73 

67.6 

0.00652 

Gravity Helmholtz {E - 20) 

3.11 

1.20 

0.0139 


For problems involving a union of a considerably large number of domains, the current method of interlacing 
all operators can be improved. In future work on fractal screens motivated by ||83l, alternative algorithms based on 
hierarchical block diagonalization via a symmetrized Schur complement ll^ may be explored specifically exploiting 
the low rank off-diagonal structure arising from coercive singular integral operators of elliptic PDOs. This is close in 
spirit to the Fast Multipole Method lIMIl . but applied to the banded representation of the singular integral operators, 
instead of discretizations arising from quadrature rules. A preliminary result in this direction is shown in the left side 
of Figure]^ 




Figure 7: Left: An illustration of a scenario that benefits from the abstraction of hierarchical matrix factorizations to our hierarchical operators. 
The Dirichlet solution of the Helmholtz equation with k = 100 and one incident source in the North Sea. Right: Idealized fluid flow around three 
obstacles. 


As illustrated in subsection 6.2 on the acoustic scattering of the Helmholtz equation with Neumann boundary 
conditions, SingularIntegralEquations. jl supports higher order diagonal singularities. Future work may ex¬ 
plore the feasibility of combining automatic differentiation and differentiation of Chebyshev interpolants to automate 
the construction of the operators with higher order singularities such that the user need only enter the fundamental 
solution with its logarithmic splitting described by ( |25] l. 

The approach developed in this article is also adaptable to other domains such as disjoint unions of circles and poly- 


25 















nomial maps of intervals and circles, see the right side ofFigurej^for an example calculated using Singular IntegralEquat ions. jl 
of idealized fluid flow over three domains: an interval, a circle and a polynomial map of an interval. To take into 
account circles, a similar analysis is straightforward with Laurent polynomials in place of weighted Chebyshev poly¬ 
nomials. However, a combined field formulation is beneficial to ensure well-conditioning when the solution of the 
exterior problem is near an eigenmode of the interior problem. Equations over maps of the unit interval and circle 
can also be reduced to numerically banded singular integral operators via approximating the map by a polynomial 
and using the spectral mapping theorem. The key formula in the Hilbert case is derived in ll68l Theorem 5.32], which 
implies that the Hilbert transform over a polynomial map of the unit interval can be reduced to a compact perturbation 
of the Hilbert transform over the unit interval. Expanding on this result, as well as adapting the procedure to log trans¬ 
forms, will be the topic of a subsequent publication. Future work may consider the use of these modified Chebyshev 
series for banded operators when two disjoint contours are in close proximity. When two or more contours coalesce, 
banded singular integral operators will depend on the ability to produce the orthogonal polynomials associated with 
that domain. Densities of the single- and double-layer potentials will have singularities on domains with cusps. Such 
an analysis is undetermined. 

As discussed in 121 , the fundamental solution of the gravity Helmholtz equation has an analogy to the Schrodinger 
equation with a linear potential. The Helmholtz equation with a parabolic refractive index shares the same analogy 
and the fundamental solution is also known mm. Parabolic refractive indices occur when considering the shielding 
of optical fibres, leading to Gaussian beams. Scattering problems in this context may shed light on the effects when 
optical fibres are occluded. Fast and accurate numerical evaluation of the fundamental solution as well as the Riemann 
function may also be possible via the trapezoidal rule. 

An important area of future research is extending the method to higher dimensional singular integral equations. 

The ultraspherical spectral method was extended to automatically solve general linear partial differential equations on 
rectangles lf84l and the ideas used to do this successfully may well translate to singular integral equations. 
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Appendix A. Proof of Theorem [19| 

To immediately satisfy the boundary conditions ( |23| l, we start with the ansatz; 

CO CO 

zo, ^o) = 1 + J] Z ~ 

1=1 j=i 

and we insert it into the integral equation (j 2 ^; 




1=1 j=i 

+ 


/' X (4 ^ ^r) ‘ - f"*' 
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dr df = 0. 


With the initial values: 
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8 i 


1 1 , 

Ao 1 —-, A 1 2 — —. A 2 2 — A, 1 /4, 

2,1 1.2 2,2 1,1/ 


and the additional values: 


(A.l) 


(A.2) 


(A.3) 


A,,i = Ai i - 0, for i > 2, A,j = 0, for i <0,j< 0. (A.4) 

the coefficients are found to satisfy in general: 

+ (4 + - ^^i-lJ -2 + = 0. (A.5) 

The growth in the constant in front of A,,y ensures that coefficients decay at least exponentially fast, hence the power 
series converges for all z and 
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To get an integral representation for the Riemann function, we start from the differential equation it satisfies after 
the change of variables u - z.- Zq and v - ^ - (q: 


d^V 

dudv 


8 i ) 


4 8i 


together with y(0,v) = V{u,0) - 1. 
Taking the Laplace transform: 


of the differential equation, we obtain: 
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Jo 


dV 1 dV i~ u\. ^ 


y(0,i) = 
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du 8i ds 

Using the method of characteristics for this first-order PDE, we obtain the general solution as: 
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The particular solution satisfying the initial condition is: 
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- {s^ - - Mij . 


(A.6) 


(A.7) 


(A.8) 


(A.9) 


(A. 10) 


Inverting the Laplace transform using the Bromwich integral, we find the solution (|150|l. 
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